﻿using System;

namespace DGEMS
{
    public static class FFT
    {

        public static double[] Calculate(double[] input)
        {
            int length;
            int bitNum;
            if (IsPowerOfTwo(input.Length))
            {
                length = input.Length;
                bitNum = Log2(length) - 1;
            }
            else
            {
                bitNum = Log2(input.Length);
                length = 1 << bitNum;
                // the items will be pad with zeros
            }

            // bit reversal
            ComplexNumber[] data = new ComplexNumber[length];
            for (int i = 0; i < input.Length; i++)
            {
                int j = ReverseBits(i, bitNum);
                data[j] = new ComplexNumber(input[i]);
            }

            // Cooley-Tukey 
            for (int fftIndex = 0; fftIndex < bitNum; fftIndex++)
            {
                int m = 1 << fftIndex;
                int N = m * 2;
                double alpha = -(2 * Math.PI / N);

                for (int k = 0; k < m; k++)
                {
                    // e^(-2*pi/N*k)
                    ComplexNumber oddPartMultiplier = new ComplexNumber(0, alpha * k).PoweredE();

                    for (int j = k; j < length; j += N)
                    {
                        ComplexNumber evenPart = data[j];
                        ComplexNumber oddPart = oddPartMultiplier * data[j + m];
                        data[j] = evenPart + oddPart;
                        data[j + m] = evenPart - oddPart;
                    }
                }
            }

            // calculate values
            double[] values = new double[length];
            for (int i = 0; i < values.Length; i++)
            {
                values[i] = data[i].AbsPower2();
            }
            return values;
        }

        private static int Log2(int n)
        {
            int i = 0;
            while (n > 0)
            {
                ++i; n >>= 1;
            }
            return i;
        }

        private static int ReverseBits(int n, int bitsCount)
        {
            int reversed = 0;
            for (int i = 0; i < bitsCount; i++)
            {
                int nextBit = n & 1;
                n >>= 1;

                reversed <<= 1;
                reversed |= nextBit;
            }
            return reversed;
        }

        private static bool IsPowerOfTwo(int n)
        {
            return n > 1 && (n & (n - 1)) == 0;
        }
    }
}
